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ABSTRACT: 

Large eddy simulation (LES) of forced, homogeneous, isotropic, two-dimensional 
(2D) turbulence in the energy transfer subrange is the subject of this paper. A dif- 
ficulty specific to this LES and its subgrid scale (SGS) representation is in that the 
energy source resides in high wave number modes excluded in simulations. There- 
fore, the SGS scheme in this case should assume the function of the energy source. 
In addition, the controversial requirements to ensure direct enstrophy transfer and 
inverse energy transfer make the conventional scheme of positive and dissipative 
eddy viscosity inapplicable to 2D turbulence. It is shown that these requirements 
can be reconciled by utilizing a two-parametric viscosity introduced by Kraichnan 
(1976) that accounts for the energy and enstrophy exchange between the resolved 
and subgrid scale modes in a way consistent with the dynamics of 2D turbulence; 
it is negative on large scales, positive on small scales and complies with the basic 
conservation laws for energy and enstrophy. Different implementations of the two- 
parametric viscosity for LES of 2D turbulence were considered. It was found that 
if kept constant, this viscosity results in unstable numerical scheme. Therefore, an- 
other scheme was advanced in which the two-parametric viscosity depends on the 
flow field. In addition, to extend simulations beyond the limits imposed by the 
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finiteness of computational domain, a large scale drag was introduced. The result- 
ing LES exhibited remarkable and fast convergence to the solution obtained in the 
preceding direct numerical simulations (DNS) by Chekhlov et al. (1994) while the 
flow parameters were in good agreement with their DNS counterparts. Also, good 
agreement with the Kolmogorov theory was found. This LES could be continued 
virtually indefinitely. Then, a simplified SGS representation was designed, referred 
to as the stabilized negative viscosity (SNV) representation, which was based on 
two algebraic terms only, negative Laplacian and positive biharmonic ones. It was 
found that the SNV scheme performed in a fashion very similar to the full equation 
and it was argued that this scheme and its derivatives should be applied for SGS 
representation in LES of quasi-2D flows. 

PACS numbers:47.25Jn 
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1. INTRODUCTION 



Homogeneous and isotropic turbulence has been a traditional idealization of real turbu- 
lent flows which are usually neither homogeneous nor isotropic. However, this idealization 
provided wealth of information on the physics of turbulence and it still remains one of the 
main tools of theoretical and numerical turbulence research (Monin and Yaglom, 1975). The 
same can be said about the dimensionality of the problem; indeed, many natural flows that 
span large number of scales possess features of both three- and two-dimensional turbulence 
and can be classified somewhere between the purely 3D and 2D extremes. Thus, even though 
the focus of this paper is isotropic 2D turbulence one should keep in mind that applications 
of the results to quasi-2D flows are thought. Quasi-2D turbulent flows are widely found in 
geophysics and engineering. Although under normal circumstances all flows are unstable to 
three-dimensional instabilities (Batchelor, 1969), there exist natural situations when a flow 
may attain a quasi-2D configuration or even become quasi-two-dimensionalized on certain 
scales. There are two major factors that may cause a flow to become quasi-2D: geometry 
of the flow boundaries and/or certain body forces (or 'extra strains') whose action leads 
to smoothing of the velocity fluctuations in a preferred direction. While in the geophysical 
context both factors are equally important (i.e., the small aspect ratio, density stratification, 
rotation), in the engineering context the second factor usually pre-dominates (for instance, 
the so called mechanism of 'magnetic friction' in magnetohydrodynamic flows with low mag- 
netic Reynolds number; Sommeria and Moreau, 1982). 

Although mathematical modeling of quasi-2D flows has important practical applications, 
particularly in the atmospheric and oceanic sciences, it has not received as much attention 
in the literature as the modeling of the 3D flows. Partly, it can be explained by the fact that 
the quasi-2D problems are less computationally intense than their 3D counterparts. Thus 
there exists a hope that in the near future, practically important quasi-2D problems can be 
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solved using DNS in which all scales are resolved (Lesieur, 1990). 

In addition, despite the specific peculiarities of quasi-2D flows related to the energy and 
vorticity dynamics, their subgrid scale representation has not received sufficient attention so 
far. There have been attempts to parameterize the SGS processes in quasi-2D flows similarly 
to those in 3D flows using Laplacian or biharmonic dissipation, the most advanced method 
being the anticipated potential vorticity method (Sadourny and Basdevant, 1985). However, 
such methods can only perform well in the vorticity dissipation subrange when energy is 
injected on relatively large scales. Being applied in the energy transfer subrange, they will 
lead to energy dissipation, contradicting basic energy and vorticity transfer dynamics of 
quasi-2D turbulence. Moreover, LES of quasi-2D flows in the energy transfer subrange has 
never been attempted, despite the fact that such flows would bear strong analogy to large 
scale oceanic and atmospheric circulation and that DNS of such flows cannot be expected in 
the foreseeable future. Thus, there exists a need to improve our understanding of the SGS 
processes in the energy transfer subrange of quasi-2D flows and to successfully simulate such 
flows when their energy sources reside in the SGS region. These both issues are addressed 
in the present paper. In the next section, the basic difficulties of the SGS representation of 
the quasi-2D flows in the energy transfer subrange are discussed. Then, the following section 
elaborates on the notion of the two-parametric viscosity and explains how this viscosity 
resolves the conflict between inverse transfer of energy and direct transfer of enstrophy. 
In section 4, advantages and deficiencies of various implementations of the two-parametric 
viscosity for LES of 2D turbulence in the energy subrange are described. In section 5, 
simplified SGS representations for LES of 2D turbulence are considered and the notion of 
the stabilized negative viscosity (SNV) is introduced. Finally, section 6 discusses the results 
and provides some conclusions. 
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2. BASIC PROBLEMS OF THE SGS REPRESENTATION OF 
QUASI-2D FLOWS IN THE ENERGY TRANSFER SUBRANGE 

Confined to two dimensions, turbulent flows become non-vortex-stretching and undergo 
dramatic structural changes (Kraichnan and Montgomery, 1980). The most profound modi- 
fications take place in the dynamics of energy and vorticity transfer. It is well known that in 
isotropic homogeneous 3D turbulence, the direct energy cascade from large to small scales 
facilitates efficient energy dissipation by molecular viscosity. This process is accompanied by 
and closely related to the production of enstrophy (mean square vorticity) through vortex 
stretching mechanism. Since in 2D flows vortex stretching can not occur, the enstrophy 
then is conserved. Thus, 2D inviscid fluids possess two nontrivial integrals of motion: the 
energy and the enstrophy. The enstrophy conservation prevents cascade of energy from 
large to small scales because such cascade would increase enstrophy (Kraichnan, 1967, 1971; 
Kraichnan and Montgomery, 1980). 

Mathematically, this important feature is illustrated by the Fj0rtoft theorem (Lesieur, 
1990). However, the direct cascade of enstrophy from large to small scales is possible. It 
results in molecular dissipation of large scale vorticity at small scales. Drawing analogy to 
eddy viscosity in 3D flows, one can infer that small scale processes in 2D flows generate 
an effective, or eddy viscosity for the vorticity of large scales. However, introducing an 
eddy viscosity concept in 2D flows seems to be intrinsically inconsistent and self-defeating 
because the dissipation of enstrophy will be accompanied by the dissipation of energy, which 
is physically incorrect. This controversy calls for modification of the eddy viscosity concept 
for quasi-2D flows; the issue in the focus of the present paper. 

More detailed consideration of the transport processes in 2D turbulence reveals that they 
depend on the wave numbers of the energy injection, kj. For k < kj-, the energy cascades up 
scales (inverse cascade), while the enstrophy flux is zero. For k > kf, the energy flux is zero, 
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but there exists the direct enstrophy flux (Kraichnan, 1971). If LES of a quasi-2D flow is 
thought, the proper SGS parameterization should depend on the wave number of the energy 
source, kf, i.e., whether kj belongs in the resolved (or explicit) or unresolved (or SGS) 
region. In the former case, a simple hyperviscous SGS representation may suffice, because 
it should only account for the enstrophy dissipation due to the direct cascade. However, if 
the forcing is located in the subgrid scales, then the hyperviscous SGS representation would 
lead to erroneous results since it will constitute energy dissipation in non-energy-dissipating 
flows. To sustain such flows, one would need to introduce a large scale energy source in the 
energy cascade subrange. A possible solution to this problem would be to replace an SGS 
forcing by a forcing located in the explicit region near k c , where k c is the cutoff wave number 
corresponding to the grid resolution. However, this solution is not only quite cumbersome 
but it also significantly distorts the explicit scales near k c . In addition, this approach is 
difficult for implementation in the physical space, particularly for bounded systems and/or 
systems with spatially non uniform energy sources. 

Another solution would be to introduce a negative eddy viscosity, first extensively dis- 
cussed by Starr (1968), as an SGS parameterization of the unresolved energy source. Recent 
studies of flows with negative viscosity were conducted by Gama et al. (1991,1994) in 2D 
and Yakhot and Sivashinsky (1985) in 3D. Although such an SGS representation could ad- 
dress the issue of the inverse energy cascade, it would not satisfy the constraint of the zero 
enstrophy flux. In addition, equations of motion with negative viscosity produce ill posed 
problems. It appears therefore that addressing the issue of SGS representation for quasi-2D 
flows in self-consistent and comprehensive way would require full consideration of energy 
and enstrophy dynamics and should be based upon the corresponding transport equations. 
Such an approach was first outlined by Kraichnan (1976) who introduced the notion of 
two-parametric viscosity. This approach and its implications will be elaborated in the next 
section. 
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3. TWO-PARAMETRIC VISCOSITY AS SGS 
REPRESENTATION OF QUASI-2D FLOWS 



Two-dimensional incompressible turbulent flows are described by the vorticity equation 

m + d( X ,y) =^C + f^ (i) 

where ( is fluid vorticity, uq is molecular viscosity and f is external force. 

The introduction of the classical eddy viscosity concept for LES with Eq. (1) implies 
that there is a distinct scale separation between the resolvable and SGS modes. Indeed, 
only if such a separation exists, the eddy viscosity would be /c-independent and a function 
of the cutoff wave number k c only. However, the assumption of scale separation fails in all 
turbulent flows, particularly in 2D flows, such that, strictly speaking, an SGS representation 
should depend on two parameters, k and k c . Such two-parametric viscosity, denoted by 
u(k\k c ), was first introduced by Kraichnan (1976). It describes the energy exchange between 
given resolved vorticity mode with the wave number k and all SGS modes with k > k c . The 
two-parametric viscosity is derived from the evolution equation for the spectral enstrophy 
density Q(k,t) = (47r) _1 /c(£(k, t)((— k, £)), where (. . .) denotes ensemble averaging: 

(d t + 2uk 2 ^)Q(k,t) = T n (k,t). (2) 
The enstrophy transfer function in (2), 7^(/c, t), is given by 



r n (M) = / ^<C(MC(q,0C(-M>^ \ . (3) 

2vr U P +q=k p z (27T) 7 



For a system in statistical steady state the two-parametric transfer T^i(k\k c ) and viscosity 
u(k\k c ) are calculated from (3) by extending integration only over all such triangles (k, p, q) 
that \k — p\ < q < k + p and p and/or q are greater than k c : 
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In a wide class of quasi-normal approximations 7^(A;|A; C ) in two dimensions is given by 



2k 



T n (k\k c ) = ^JJ ®- klP ,q(p 2 ~ <? 2 ) sin a 



2 2 



2 2 2 2 

; ' ~ q -n(q)n(k) + k ~l n( P )n(k) 



dpdq, (5) 



A; 3 g 3 A; 3 p 3 

where Q^k,p,q i s the triad relaxation time. Also, the angle a is formed by the vectors p and 
q, and f denotes integration over the area defined above. 

Different spectral closure models provide different specification of Q-k,p,q- Chekhlov 
et al. (1994) compared u(k\k c ) calculated from their DNS data with 512 2 resolution with 
those evaluated by Kraichnan (1976) using his Test Field Model (TFM) and obtained from 
the Renormalization Group (RG) theory (see Appendix). The results of this comparison 
are shown in Figs. 1 and 2. Figure 1 presents the DNS-inferred normalized two-parametric 
viscosity, 

N(k/k c ) = v(k\k c )/\v(0\kc)\, (6) 

along with the TFM- and RG-based analytical predictions. The results are in a very good 
agreement with each other over the entire energy transfer range, up to the wave numbers 
close to k c , where the DNS data saturates, while TFM and RG curves exhibit sharp cusp. 
This theoretical cusp is due to the fact that as k — > k c , more and more elongated triads with 
either p or q ^ k c become involved in the energy exchange between the mode k and the 
subgrid scale modes. The contribution from these triads to the energy exchange near k c is 
very significant and results in the cusp behavior. However, in finite box DNS with large-scale 
energy removal, the energy of small wave number modes is reduced and u{k\k c ) is expected 
to saturate near k c . Indeed, when the RG-based v(k\k c ) was re-calculated based upon the 
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DNS energy spectrum, the unnormalized DNS- and RG-based two-parametric viscosities 
were found to be in very good agreement for all wave numbers, as shown in Fig. 2. 

Figures 1 and 2 show that for the large scale modes, for which k <tik c and scale separation 
exists, the effect of the SGS modes is represented by a negative and constant viscosity, such 
that these modes gain energy from their SGS counterparts by means of the inverse transfer. 
On the other hand, u(k\k c ) > for k — > k c such that the modes close to k c lose their energy 
to the SGS modes. The difference between the large scale gain and the small scale loss is 
equal to e, the rate of the energy input due to the forcing / in Eq. (1); see also Eq. (13) 
below. If enstrophy balance is considered, recall that the enstrophy transfer is most efficient 
at small scales, such that the resulting balance for the resolvable scales turns out to be zero 
[Kraichnan, 1971; see also Eq. (18) below], i.e., the enstrophy is conserved. This explains 
how the two-parametric viscosity resolves the controversy of the inverse cascade of energy and 
conservation of enstrophy in the energy subrange of 2D flows. It appears therefore that the 
only physically correct way to represent SGS processes in 2D turbulence would be through 
the two-parametric viscosity. Such an approach has become quite popular in simulations of 
3D flows (Domaradzki et al., 1987; Galperin and Orszag, 1993) but has not yet been fully 
explored for LES of quasi-2D flows. The implementation of the two-parametric viscosity for 
LES of 2D turbulent flows, the arising problems, their solutions and results are described in 
the following sections. 
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4. IMPLEMENTATION OF THE TWO-PARAMETRIC 
VISCOSITY FOR LES OF 2D TURBULENCE 



To test the two-parametric viscosity-based SGS parameterization in the energy transfer 
subrange, a series of LES of 2D turbulence in Fourier space was designed. These LES were 
based upon Eq. (1) in which all subgrid scale processes including the forcing were represented 
by the two-parametric viscosity u(k\k c ), 

+ / ^c( P )C(k - p)t^2 = K*i*c)* 2 c(k), o < k < k c . ( 7 ) 

ot J\p\,\k-p\<k c P z (2?rr 
It is important to reiterate that in LES of 2D turbulence in the energy subrange, the source 
of energy resides on the unresolved scales, such that Eq. (7) appears unforced. However, 
as was explained earlier, the negative part of u{k\k c ) serves as the only energy source for 
the resolved modes. In the course of the present LES it was found that numerical results 
critically depend on the way u{k\k c ) is introduced in the solver. Thus, a series of simulations 
was designed with the purposes of understanding the nature of the problems associated with 
the implementation of the two-parametric viscosity and of identifying the most viable and 
robust ways to use this viscosity in LES of 2D flows. 

4.1. Description of numerical method 

The numerical solver used in the present calculations was based upon Fourier- Galerkin 
pseudo-spectral formulation (Orszag, 1969) in the periodic box [0, 2tt] x [0, 2n] and was 
essentially the same as that utilized in the preceding DNS (Chekhlov et al., 1994). The 
present LES employed 162 2 resolution with complete dealiasing based on the 2/3-rule; the 
cutoff wave number was set at k c = 50, which is about half of the resolution used in DNS 
of Chekhlov et al. (1994). The preprocessing was done using the second order Runge- 
Kutta scheme, while marching in time was accomplished using an implicit, second order, 
stiffly stable Adams scheme (Karniadakis et al., 1991). The initial flow field was set to zero 
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everywhere except for a narrow band of wave numbers in the middle part of the spectrum 
where it was assigned random Gaussian distribution. The DNS inferred value of e was about 
5.19 x 1CP 10 ; the same e was used in LES. The time step in LES was set to St = 0.5 satisfying 
both convective and viscous necessary conditions for linear stability. Based upon the size 
of the largest energy containing eddy with the wave number fc m i n , 27r//c m j n , and the total 
energy of the steady-state E(t), the maximum large scale eddy turnover time defined as 
Ttu = 27r/(/c m j n V / 2£') was about rt u ~ 3600 for LES of cases 1 and 2 below, where k m [ n = 1, 
and Tt u — 900 for cases 3,4, and 5, for which k m [ n = 4. 

4.2. Case 1. Flow independent v{k \ k r ) 

With u{k\k c ) known and flow independent, Eq. (7) can be solved directly. According 
to Fig. 1, u(k\k c ) can be obtained from DNS or from some statistical theory of turbulence. 
Thus, in the first LES numerical solver for Eq. (7) utilized u(k\k c ) derived from the renor- 
malization group (RG) theory of turbulence (see Appendix for the details). As shown in the 
Appendix, 

u(k\k c ) = 0.327e 1 / 3 A; £ : 4/3 iV(A;/A; c ), (8) 

where N{k/k c ) is given by (6), and since e is a constant, v(k\k c ) is a function of k and 
k c only. It was assumed that thus defined u{k\k c ) would be capable of supporting inverse 
energy cascade with constant SGS energy input e. To verify this assumption, one needs to 
examine the evolution of total energy and enstrophy of the resolved modes, E(t) and Q(t), 
respectively. By definition, E(t) = Jq c E(k,t)dk, where E(k,t) = (4vrA;)" 1 (C(k, i)C(-k, t)) 
is the spectral energy density, and Q(t) = Jq c Q(k,t)dk. The basic requirement to LES 
would be that E(t) and Q(t) of LES have the same behavior as those derived from Eq. (1) 
for which the evolution laws are E(t) oc et and Q(t) = const, due to the conservation of the 
inviscid integrals for E(t) and Q(t) (recall that in the energy subrange of 2D turbulence the 
rate of the enstrophy flux rj = 0). Figures 3a,b show that in the first LES, both E(t) and 
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Q(t) exhibit nonlinear growth indicating that not only the rate of the energy transfer to the 
resolvable scales e LE s 7^ const, but also the rate of the enstrophy transfer % ES 7^ 0. Figure 
4 shows that instantaneous spectrum E(k,t) also reveals tendency to growing up with time 
without stabilizing around any universal distribution. At some point of its evolution, E(k, t) 
crossed the Kolmogorov law 



where Cj( m 6 is the Kolmogorov constant, and e was close to its prescribed DNS value. 
However, at larger times the Kolmogorov scaling (9) was lost, while €les kept growing. The 
roots of the problem are revealed when one calculates the rate of the energy input into all 
resolved modes, ^les- The energy equation derived from the definition of spectral energy 
density and Eq. (7) yields 



Since in the first series of LES, v{k\k c ) depends on k and k c only, and E(k,t) is a dynamic 
variable that depends on the evolution of the flow field, cles(^) a l so turns out to be time 
dependent. This is not only in direct conflict with the requirement that €les = e =const, but 
also leads to a positive feedback between the energy input and total energy of the system, 
which results in numerical instability. To correct this problem, one must ensure that 6les = 
const. This can be achieved by allowing u(k\k c ) to become time dependent and related to 
resolved variables. The philosophy of using actual flow field characteristics to determine 
energy input and dissipation would be analogous to a standard practice of 3D LES. 

Inspection of Eq. (10) reveals that cles(^) is proportional to the total enstrophy of 
the resolvable field, such that stipulating 6les = const would require u(k\k c ) to become time 
dependent and inversely proportional to Q(t) (the time dependency of u(k\k c ) will be implied 
in the following discussion but suppressed in notations). To find the explicit dependency of 



E(k) = C K e 2 / 3 k- 5 / 3 , 



(9) 




k, 



(10) 
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u(k\k c ) on Q(t), let us assume that u{k\k c ) can be represented by 

v(k\k c ) = F(H)N(k/k c ), (11) 

where F(-) is some function of Q(t) and N(k/k c ) is still given by (6). According to Fig. 1, 
N(k/k c ) can be split into negative and positive terms, 

N(k/k c ) = -1 + <p(k/k c ), <p(k/k c ) > 0, < k/ k c <l. (12) 

One can now calculate integral in (10) using Eqs. (11) and (12): 

e LE s = 2F(n) [ ° E(k,t)k 2 dk-2F(Q) ! ° <f)(k/k c )E(k,t)k 2 dk 
Jo Jo 

= 2F(U)U(t) - 2F(U) / " <p(k/k c )E(k,t)k 2 dk. (13) 

JO 

Integration of (13) for the Kolmogorov spectrum with (p(k/k c ) evaluated from the RG theory 
(see Appendix) yields 

e LES ^ O.SF(U)U(t). (14) 
Equation (14) shows that to satisfy the requirement €les = e = const, one has to impose 

F(H)=e/[0m(t)], (15) 

such that the two-parametric viscosity (11) becomes 

v(k\kc) = - N(k/k c ) = L— + — ^Mk/k c ). (16) 

The first term in the right hand side of (16) accounts for the SGS energy input while the 
second term represents the high wave number dissipation as k — > k c . As was argued earlier, 
to ensure €les = const, the energy source term must be time dependent and inversely 
proportional to Q(t). Thus, the negative feedback between energy input and enstrophy 
of the resolved modes is the mechanism that stabilizes numerical process. Note that the 
SGS formulation based upon Eq. (16) complicates Eq. (7) because its right hand side now 
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depends on the functional of the solution, Q(t). However, on the one hand, it is clear from 
the presented analysis that LES of 2D turbulence based upon Eq. (7) is impossible if v{k\k c ) 
depends on k and k c only. On the other, the SGS representation (16), though complicated, 
is in line with the eddy viscosity approach, in which the eddy viscosity coefficient is usually 
solution dependent (Smagorinsky, 1963,1993; Yakhot and Orszag, 1986; Yakhot et al., 1989). 

4.3. Case 2. Flow dependent v(k \ k r ) with no large scale drag 

The SGS formulation (16) was used in the second LES and considerable improvement 
over Case 1 was observed. Figures 5a,b show that up to the simulated time t = 3r^ u , 
E(t) grows linearly, while Q(t) attains a constant value. Figure 6 shows that during the 
same t, E(k,t) quickly approaches steady state Kolmogorov distribution (9). However, for 
t > 3rt u the flow field undergoes irreversible modifications; the behavior of E(t) and Q(t) 
changes, while E(k,t) begins to deviate from the Kolmogorov law. All these changes reflect 
the basic problem of the present LES that simulates the behavior of an infinite system in 
a finite computational box (Smith and Yakhot, 1993a,b). In this box, the smallest wave 
number modes become energy saturated at t ~ 3r^ M , and, if the energy of these modes is not 
removed, they begin to alter the behavior of the entire flow field. Therefore, to extend LES 
beyond t ~ 3r^ u , one needs to prevent the accumulation of energy at the lowest modes, which 
was accomplished in LES of Case 3. Note however that by the time t ~ 3r^ u the inverse 
cascade swept through all the resolved modes such that they became energy saturated and 
attained the steady state. Therefore, one should expect that in LES with t > 3r^ M both E(t) 
and Q(t) remain nearly constant. 

4.4. Case 3. Flow dependent u(k \ k r ) with large scale drag 

The simplest way to withdraw energy from the lowest modes would be by mere setting to 
zero the amplitudes of those modes. However, such a "chopping" alone is known to produce 
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unsatisfactory results (Browning and Kreiss, 1989). Therefore, in addition to the chopping, 
one needs to introduce a mechanism that would account for the energy exchange between 
the resolved modes and the low wave number modes excluded in LES. Such a mechanism, 
the large scale drag, was introduced in this study in analogy to the two-parametric viscosity. 
However, the derivation of this large scale drag is beyond the scope of this paper and will 
be presented elsewhere. 

The large scale drag was implemented in the third LES, whereas the amplitudes of all 
modes with k < k m [ n = 4 were set to zero. As was explained in section 3 and shown in 
Fig. 2, such a chopping results in flattening of the cusp in N(k/k c ) as k — > k c , such that 
this function had to be recalculated which in turn led to modification of the coefficient in 
(14) from 0.8 to 0.87. The results of the third LES are shown in Figs. 7a,b and 8a,b. 
One can see that this simulation could virtually be extended indefinitely, with E(t) and 
Q(t) slightly oscillating around their steady state values (the source of these oscillations is 
probably the self- adjustment of the numerical scheme to the mismatch between the small 
scale energy forcing and large scale withdrawal). The instantaneous energy spectrum, Fig. 
8a, exhibits steady and nearly perfect Kolmogorov scaling. Since the large scale drag enables 
one to dramatically increase the integration time in LES, it will be retained in all further 
simulations. Note however that these simulations will pertain to steady state rather than 
time developing flows. 

An important characteristic of 2D turbulence in the energy transfer subrange is the 
energy flux Ti E (k) = Jq T^r^nT^dn which theoretically should be equal to e for any k. Note 
that the enstrophy transfer function Tq(ti) is defined by (3) where integration is extended 
over all triangles k + p + q = including the SGS ones for which k and/or p and/or 
q > k c . However, if n^(A;) is calculated for LES results, then the SGS triangles should be 
excluded such that the integration area in (3) is reduced. On the other hand, the energy 
flux into the interval < n < k in LES consists of two contributions, Jq T^(n)n^ 2 dn 
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(where 7^ implies the reduced integration area in (3)), and the direct flux from the SGS 
modes, —2 Jq u(n\k c )n 2 E{n)dn. Figure 8b shows that thus calculated remains nearly 

constant for k > 16 and is equal to about 5 x 1CP 10 , which is very close to the specified 
value of e obtained from DNS by Chekhlov et al. (1994). The decrease of Ti^(k) at small k 
is due to the large scale energy removal. 

4.5. Case 4. Flow dependent energy input with flow independent dissipation 

It would be tempting to simplify (16) by relaxing the time dependency in the dissipation 
term. It is not clear a priori whether or not this time dependency is critical, and to find 
out about it a fourth LES was conceived in which SGS representation (16) was modified by 
replacing e/0.8 VL{t) in the dissipation term by the RG derived expression 0.327e 1 / 3 A; c 4 ^ 3 . 
Figures 9a,b and 10a,b,c show that this simplified SGS scheme performs in a very robust 
way with no oscillations at all for a relatively long time, t ~ 30t£ W . Similarly to Case 3, 
there exists a mismatch between the small scale forcing and large scale energy removal, but 
since the dissipation in Case 4 cannot self-adjust, the solution begins to deteriorate when 
this mismatch accumulates significantly. Still, Case 4 LES could be extended to many more 
turnover times than the corresponding DNS, the result quite remarkable for its own sake. 
During this time, the instantaneous energy spectrum, Fig. 10, exhibits steady and nearly 
perfect Kolmogorov scaling almost indistinguishable from that in Fig. 8. Figure 10b shows 
that the Kolmogorov constant Cx in (9) is about 5, in good agreement with the RG derived 
value of 5.12 (see Appendix). As in Case 3, the energy flux He(^) shown in Fig. 10c is 
almost constant for k > 16 and is about —5 x 10~ 10 . 

Although the SGS formulation (16) is relatively easy to implement in spectral LES, 
in practical applications, particularly in the physical space, it would be far more useful to 
approximate N(k/k c ) analytically. An additional benefit of such an approach would be a 
possibility to carry out further analytical studies of this SGS representation; the direction 
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pursued in the next section. 
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5. STABILIZED NEGATIVE VISCOSITY (SNV) FORMULATION 



For practical implementation of SGS formulation (16) it is convenient to approximate 
N(k/k c ) in (6) by a series in powers of k 2 . It was found that even the first two terms of this 
series, 

N(k/k c ) ~ -1 + a(k/k c f, (17) 

where a is a constant, are sufficient to perform successful LES of 2D turbulence. To find a 
recall that representation (17) must ensure zero enstrophy transfer in the energy subrange, 

r/ LES = 2 [ C u{k\k c )E{k, t)k A dk = 2F(Q) [ ° N{k/k c )E{k, t)k 4 dk = 0. (18) 

Jo Jo 

Substituting (17) into (18) and assuming that E(k, t) is Kolmogorovian, one finds that a — | 
such that (13) becomes 

e LES = -2F(H) f kc [-l + l(k/k c f]E(k,t)k 2 dk 
Jo 5 

= 2F(H)H(t) - —F(U)k- 2 P(t), (19) 

where P(t) = Jq c E(k,t)k 4: dk is the total palinstrophy of the resolved modes. For Kol- 
mogorovian E(k,t), (19) can be integrated analytically to yield P(t) = ^^l(t)k 2 and 

m - ¥s m (20) 

Equation (20) completes the two-term SGS parameterization for LES of 2D turbulence in 
which the right hand side of Eq. (7) takes the form 



2,— 25 e 



v(k\k c )k z ((k) = -=- 



1 + sU 



k- 2 ' 



fc 2 C(k). (21) 



18 Q(t) 

As in (16), the first term in the right hand side of (21) accounts for the energy flux from the 
unresolved modes and the inverse cascade of this energy, while the second term represents 
the energy dissipation near the cutoff. Using (21) and following the philosophy of Cases 
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3 and 4 LES, two more simulations were designed, with flow dependent and independent 
dissipation term in (21). 

5.1. Case 5. Two-term LES with flow dependent dissipation 

Simulations performed with the formulation (21) slightly adjusted to account for the 
finiteness of the computational domain are shown in Figs. lla,b and 12. They exhibit very 
little difference compared to the third LES that employed the full curve N(k/k c ) given by 
(6) and shown in Figs. 7a,b and 8. One infers therefore that (21) is a viable two term SGS 
representation for LES of 2D turbulence; obviously, (21) is significantly simpler than (16). 

5.2. Case 6. Two-term LES with flow independent dissipation 

For practical purposes, it would be most appealing to use formulation (21) with the 
dissipation term constant. A numerical experiment analogous to that of Case 4 LES was 
conducted with the energy source in (21) not changed but in the dissipation term, Vt{t) was 
replaced by its value calculated for the Kolmogorov spectrum (9). Such an approach yields 
the dissipation term in (21) in the form 

Ak 4 ((k), (22) 

where A is a constant given by 

A = O.Slle 1 / 3 ^ 1073 , (23) 

which corresponds to (7^=5.8. The results of this case 5 LES are presented in Figs. 13a,b 
and 14a,b; there is a very good agreement with the corresponding results obtained with the 
full curve N(k/k c ) up to t ~ 40t£ U after which, similarly to Case 4, the solution begins 
to deteriorate. The Kolmogorov constant shown in Figs. 14b is close to 5.8, which is 
consistent with derivation of (23). One can therefore infer that (21), (23) is probably the 
simplest SGS representation for LES of 2D turbulence possible. 
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Further advantages of the SGS representation (23) are revealed when LES of quasi-2D 
turbulence is sought in the physical space, where Eq. (21) combined with (23) leads to the 
following SGS representation: 

-IX A ^ ] )- A w^ (24) 

1 J 

where 

A ^lwy (25) 

A A = A = 0.511e 1 / 3 (A/27r) 10 / 3 = const, (26) 

and where fi(x) denotes the enstrophy averaged over an area adjacent to the grid cell, and 
A is the grid resolution (note that the Laplacian term in (24) is written in the conserva- 
tive form). Equation (24) thus includes two terms, the negative Laplacian and positive (in 
the sense of dissipation) biharmonic, and structurally resembles the Kuramoto-Sivashinsky 
equation widely known from combustion theory (Sivashinsky, 1977) and flows with chemical 
reactions (Kuramoto and Tsuzuki, 1976; Kuramoto, 1984). However, the SGS representa- 
tion (24)- (26) combined with the explicit equation for the resolved scales produces far more 
complicated equation than the Kuramoto-Sivashinsky equation because generally its coef- 
ficients are not constant but, as in the eddy viscosity approach, are functions of the flow. 
There have been previous attempts to use formulation similar to (24)- (26) but with constant 
coefficients (dubbed the Kuramoto-Sivashinsky-Navier-Stokes equation, see Gama et al., 
1991) to perform LES of 2D turbulence. However, they were not overly successful even in 
reproducing the Kolmogorov spectrum, mostly because they used constant "eddy viscosity" 
coefficients. Since, on the one hand, SGS representation (24)-(26) includes a negative Lapla- 
cian viscosity term and a positive, stabilizing, dissipation term, but, on the other hand, it 
is quite different from the Kuramoto-Sivashinsky equation in which "viscosity" coefficients 
are constant, it will be referred to as the Stabilized Negative Viscosity (SNV) formulation. 
The SNV formulation is expected to be particularly useful in simulations of atmospheric 
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and oceanic flows where large scale motions are quasi-two-dimensional by their nature while 
small scale forcing is significant (Monin, 1986). 
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6. DISCUSSION AND CONCLUSIONS 



LES of homogeneous and isotropic 2D turbulence in the energy transfer subrange and 
appropriate SGS representation have been the central subject of this paper. Conservation 
of energy and enstrophy in 2D turbulent flows leads to coexistence of two spectral transfer 
processes, upscale for energy and downscale for enstrophy. Conventional eddy viscosity 
formulations are purely dissipative and fail to accommodate both transfers simultaneously. 
It is argued that proper SGS parameterization for LES of 2D flows is given by the two- 
parametric viscosity u(k\k c ) introduced by Kraichnan (1976) that accounts for the energy (or 
enstrophy) exchange between given resolved and all SGS modes and which includes negative 
and positive branches; the negative branch of u(k\k c ) represents the unresolved, small scale 
forcing and inverse cascade of energy, while the positive one represents the dissipation. It was 
shown that the negative and positive parts of the two-parametric viscosity play vital role in 
ensuring that all conservation and spectral transfer laws of 2D turbulence are satisfied. Then, 
v(k\k c ) was used in LES of 2D turbulence in the energy transfer subrange where the negative 
part of v{k\k c ) was the only energy source. The sensitivity of numerical results to the way 
of implementation of v(k\k c ) in numerical schemes was studied. It was found that if v{k\k c ) 
is specified as a flow independent parameter then a positive feedback is established between 
the forcing and the total energy of the system leading to numerical instability. Thus, another 
scheme was designed in which u(k\k c ) was flow dependent but the rate of the energy input 
was kept constant. This LES exhibited a very stable behavior consistent with analytical 
theories and DNS; Kolmogorov scaling was evident and robust, and all conservation and 
spectral transfer laws were fulfilled. Then, a simplified SGS representation was advanced, in 
which only two terms were retained: one negative and proportional to A; 2 , and one positive 
and proportional to A; 4 . The viscosity coefficient in the negative term served as the energy 
source and had to be flow dependent, to ensure that the energy input remains constant. 
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However, since the A; 4 term is mostly active at high wave number region near the cutoff wave 
number k c , and its main role is energy and enstrophy dissipation at k — > k c , this term did 
not need to be flow dependent and could essentially be kept constant and evaluated from 
analytical theories of 2D turbulence. Indeed, LES with constant and flow dependent A; 4 
terms gave very similar results. 

The application of the two term parameterization to simulations in physical space results 
in the so called Stabilized Negative Viscosity (SNV) representation which includes negative 
Laplacian and positive biharmonic viscosities. Numerical implementation of the scheme 
with the constant biharmonic viscosity is obviously much simpler than that with the flow 
dependent viscosity. 

On the one hand, the negative viscosity term is essential in SNV scheme; on the other, 
this scheme substantially differs from the Kuramoto-Sivashinsky equation with its constant 
viscosity because in SNV, the negative viscosity not only is not constant but is a nonlinear 
functional of the solution. In fact, the SNV representation is a peculiar case of the eddy 
viscosity approach; one normally expects that if all scales of a turbulent flow are resolved then 
eddy viscosity would become equal to the molecular viscosity which is true in 3D turbulence. 
In 2D turbulence, however, resolution of all scales must be accompanied by restoration of the 
explicit small scale forcing which would result in disappearance of both negative Laplacian 
and positive biharmonic viscosities and appearance of one positive, Laplacian, molecular 
viscosity. 

It is believed that the SNV representation should be especially useful in quasi-2D flows 
in which considerable amount of energy resides on the unresolved scales. Such flows are 
typical in geophysical fluid dynamics where, due to the topographic and other constraints, 
flows are quasi-2D and where the small scale forcing is the predominant source of energy. 
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8. APPENDIX: DERIVATION OF FLOW 
PARAMETERS BASED ON RG THEORY 



Fourier transformed Eq. (1) reads 

ah) = G°(k) I k g x 2 qC( ^ )3 ~' ) ^ + G°m (ai) 

where k = (k, cu), q = (q, Q), and G°{k) = (iu> + z/qA: 2 )" 1 is the bare propagator. The 
random, zero- mean, white in time Gaussian stirring force / in the right hand side of (Al) 
accounts for the steady energy input localized in the vicinity of a given wave number kj. In 
the RG formalism (Yakhot and Orszag, 1986; Staroselsky and Sukoriansky, 1993) the effect 
of a localized random stirring on the bare equation of motion is equivalent to the effect of a 
spatially-distributed forcing with the correlation function 

< /(k, w)/(k', J) >= D{k)(2nf5(k + k')5(u + J) (A2) 

on the renormalized equation. Here, D{k) = 2Dok~ y+2 , where y and Dq are different for 
k <^ kj- and k ^> kj-. The region k <C kj corresponds to the inverse cascade of energy for 
which y = 2 and Dq oc e, e being the constant energy injection rate; the region k 3> kf 
corresponds to the direct cascade of enstrophy with the constant rate 77; in that region 
Dq oc i] and y = 4. Note that in the framework of the 2D RG theory, the regions k kf 
and k ^ kj- are treated as separate asymptotic regimes. Equation (Al) is defined in the 
interval < k < Aq. 

Fourier coefficients of the vorticity field can be separated into small scale, or fast modes, 
C > (^), for which Aq — SAq < k < Aq, and large scale, or slow modes, C < (^) ) f° r which 
< k < Aq — 5Aq, respectively. The RG procedure detailed in Forster et al. (1977) and 
Yakhot and Orszag (1986) consists of successive elimination of infinitesimal shells of the 
small scale modes from the equation of motion for the large scale modes. Essentially, ( field 
decomposed into C > (A ; ) an d C < (A ; ) parts is substituted into Eq. (Al); then, the one-loop 
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approximation is invoked and the average is taken over the small scale modes. Equation (Al) 
is then rewritten in terms of C < (^) and is defined in the abridged interval < k < Aq — 5Aq; 
it now contains a new term represented by the correction to the Green function, 

= /> fcV-(k. q) 2 Q,--2_A Git-e)MQl>D«->»$ s , (A3) 

that accounts for the effect of the unresolved scales. Here, J > denotes an integration over 
the band of wave numbers being removed. 

The integral in (A3) is calculated in the limit k — > 0, uo — > 0. The 0(k 2 ) terms provide 
correction to uq. Iterating the scale elimination procedure, one can remove a finite band of 
modes and obtain a differential equation for the parameter of renormalization z/(/c), 

dv D 



dk 167iu 2 k e+1, 

where e = 4 + y — d. For the energy inertial subrange considered below, y = 2 and e = 4. 
The solution to Eq. (A4) is 

3D 



(A4) 



U(k) = UQ < 1 + 



1/3 

(A5) 



Assuming that fc/Ag 1 and the molecular viscosity can be neglected [uq <^ v(k)], one 
finds: 

"(*)= f^v /3 ^ 4/3 - (A6) 



v 64tt^ 

It is argued (Forster et al., 1977; Yakhot and Orszag, 1986) that iterated indefinitely, the 
RG procedure of small scale elimination converges to a fixed point solution whereas ((k) is 
described by the Langevin equation 

C(k,cj) = G(k,cj)f(k,cj), (A7) 

where G = [icu + //(A;)^ 2 ]" 1 is the renormalized propagator. The existence of the fixed point 
for RG procedure in 2D has been proven by Staroselsky et al. (1995) for e I. Although the 
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feasibility of continuation of these results into the region of large e is still an open question, 
it is assumed here that the fixed point solution (A7) exists at e = 4. 

Equation (A7) allows one to calculate the vorticity correlator, U (k, uo) = (C(k, t)((— k, t)), 
as well as the kinetic energy spectrum, 



Although u{k) describes renormalization of the bare viscosity uq, it is not what is often 
comprehended as an eddy viscosity. This is merely a response function of the nonlinear 
dynamical system described by the renormalized Eq. (A7). As such, it allows for calculating 
vorticity correlator and energy spectrum but does not directly relate to enstrophy and energy 
transfer and dissipation. Furthermore, [z/(/c)/c 2 ] -1 can be viewed as a characteristic time scale 
of information loss at given k caused by nonlinear scrambling of all other modes (Dannevik et 
al., 1987). Therefore, z/(/c) has a meaning of an eddy damping parameter and is substantially 
one-point turbulence characteristic. 

To analyze energy and enstrophy transfer, one needs to consider two-point characteristics 
that account for interaction between a given explicit mode k < k c and all modes k > k c ; k c 
is identified with the moving dissipation cutoff. For this purpose, following Dannevik et al. 
(1987) the energy evolution equation should be derived at lowest nontrivial order of nonlinear 
coupling using fully renormalized propagator G(k). The resulting dynamical closure is similar 
to the Eddy-Damped, Quasi-Normal, Markovian (EDQNM) approximation where the lowest 
order RG analysis is invoked to obtain the eddy damping function. The energy equation 
then reads 




(A8) 




(A9) 



where 
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1,2 2 t.2 2 

' q -E(q, t)E(k, t) + —^E(p, t)E(k, t) 



sin a, (A10) 



kq kp 

and where vj. = u(k)k 2 , a is an angle opposite to the vector k in the triangle k + p + q = 0, 
and the integration domain A is defined by the triangular inequalities \k — p\ < q < k + p. 
The function u{k) given by (A6) is used to compute the relaxation time 0^p q {t) = (1 — 

Following Kraichnan (1976) one can now define the two-parametric, or effective eddy 
viscosity at wave number k in terms of the energy transfer from all subgrid scale modes with 
k > k c to the given explicit mode k, 

v{k\k c ) = -T(k\k c )/[2k 2 E(k)], (All) 

where 

T(k\k c ) = J J^T(k,p,q)dpdq, k < k c , (A12) 

and where the integration is extended over all p and k such that p and/or q > k c . Assuming 
that the limit t — > oo in (A10) is considered the time argument in T(k,p,q) has been omitted. 
In that limit, 6^ = (v^ + v p + 

For k <^ k c the two-parametric viscosity u{k\k c ) can be calculated analytically (Kraich- 
nan, 1976). In this case, the triangular inequality becomes \p — q\ < k q. Therefore, all the 
quantities that enter T(k,p, q) can be expanded in powers of p — q. Then, the p integration 
can be performed resulting in 

1 f°° d 

v(k\k c ) = -J^ 9 kqq — [qE(q)}dq, k < k c . (A13) 

Substitution of (A6) and (A8) into (A13) gives the asymptotic eddy viscosity for the 
largest scales, 

"(«=-KS) 1/3 ^ 4/3 (Ai4) 
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For arbitrary k, u{k\k c ) was calculated via numerical integration of (A10); the resulting 
normalized two-parametric viscosity N(k/k c ), Eq. (6), is shown in Fig. 1. This function is 
negative for k <C k c and positive for k — > k c . 

Noting that in the energy transfer subrange the energy injection rate e is equal to the 
rate of energy transfer from all the subgrid modes k > k c to all explicit modes k < k c , one 
can find the relation between the forcing amplitude Dq in (A2) and e. Indeed, as follows 



Substituting (A8) and (A14) into (A15) and performing numerical integration one finds 



Using (A16) one can now calculate the Kolmogorov constant Cj{ in (A8) and the numerical 
factor in (A 14), which are respectively, 



from (All) and (A12), 




(A15) 



D ~ 63e. 



(A16) 



E(k) = CVe 2/3 AT 5/3 , C K ~ 5.12, 



(A17) 



and 



v 



(0|fc c ) = -0.327e 1 / 3 A; c : 4/3 . 



(A18) 



These results have been used in Eqs. (8) and (14). 
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9. LIST OF THE FIGURE CAPTIONS 



Figure 1. Normalized two-parametric eddy viscosity from DNS (dots), from TFM (dashed line), 
and from RG (solid line) (from Chekhlov et al., 1994). 

Figure 2. Actual two-parametric eddy viscosity from DNS (dots) and from RG (solid line). In RG 
calculations, the energy spectrum for k < 5 was corrected in accordance with the DNS 
data (from Chekhlov et al, 1994). 

Figure 3. The evolution of the total energy E(t) (a) and total enstrophy Q(t) in Case 1 LES. Figure 
3a also shows the evolution of E(t) with the energy of the 1st, 2nd, 3rd, 4th, 5th, 6th 
and 7th modes removed. 

Figure 4. The evolution of the instantaneous energy spectrum for t/rt u = 0.56, 1.11, 1.67and2.78. 
The solid line shows the Kolmogorov -5/3 slope. 

Figure 5. Same as Fig. 3 but for Case 2 LES. 

Figure 6. Same as Fig. 4 but for Case 2 LES. Note that after t/rt u ^ 2 all instantaneous profiles 
E(k,t) become close to Kolmogorov law (9). 

Figure 7. Same as Fig. 3 but for Case 3 LES. Because the amplitudes of the first four modes are 
set to zero, only the evolution of E(t) with the energy of the 4th, 5th, 6th and 7th modes 
removed is shown. 

Figure 8a. Same as Fig. 6 but for Case 3 LES. Note that since the Kolmogorov scaling is attained 
after about t/rt u = 2, only time average energy spectrum is shown. 

Figure 8b. The energy flux U E (k) for Case 3 LES. 

Figure 9. Same as Fig. 7 but for Case 4 LES. 

Figure 10a. Same as Fig. 8 but for Case 4 LES. 

Figure 10b. The time averaged Kolmogorov constant Cx for Case 4 LES. 
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Figure 10c. The energy flux Il E (k) for Case 4 LES. 

Figure 11. Same as Fig. 7 but for Case 5 LES. 

Figure 12. Same as Fig. 8 but for Case 5 LES. 

Figure 13. Same as Fig. 7 but for Case 6 LES. 
Figure 14a. Same as Fig. 8 but for Case 6 LES. 

Figure 14b. The time averaged Kolmogorov constant Cj( for Case 6 LES. 
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